function y = impedance(m_max, CONSTS, plot_data)

    a = CONSTS.a;
    k0 = CONSTS.k0;
    d = CONSTS.d;
    c = CONSTS.c;
    eps = CONSTS.eps;
    eta = CONSTS.eta;
    eps_a = CONSTS.eps_a;
    Z0 = 4*pi/c;
    phi_vec = 0;    
    coef_sgs_to_si = 9e+11;
    
    % input impedance obtained numerically
    I_sum_0 = func_I_sum(phi_vec, m_max, CONSTS); %func_a0(CONSTS);
    y = coef_sgs_to_si/I_sum_0;
    
    % input impedance in the case of a homogeneous magnetoplasma
    h = k0*(abs(eps*eta))^(1/4)*((1-1i)/sqrt(2)); 
    I_sum_odn = -((1i*pi*h)/(Z0*k0*log(4*a/d)))*((cos((pi-abs(phi_vec)).*h.*a))./(sin(pi*h*a)));
    y1 = coef_sgs_to_si/I_sum_odn;
    
    % input impedance obtained analytically
    eps_ef = ((2i*eps_a)/((1+(eps_a)^2/(abs(eps*eta)))*(sqrt(abs(eps*eta))))) - 4;
    h1 = k0*(abs(eps*eta))^(1/4)*((1+1i)/sqrt(2*eps_ef));
    I_sum_analytical = -((1i*pi*h1)/(Z0*k0*log(4*a/d)))*((cos((pi-abs(phi_vec)).*h1.*a))./(sin(pi*h1*a)));
    y2 = coef_sgs_to_si/I_sum_analytical;

    if(plot_data)
        fprintf('Z = %f + i%f\n', real(y), imag(y));
        fprintf('Z_an = %f + i%f\n', real(y2), imag(y2));
        fprintf('Z_odn = %f + i%f\n', real(y1), imag(y1));
    end

end